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SUMMARY 


The fact that people think or behave differently from 
one another is rooted in individual differences in 
brain anatomy and connectivity. Here, we used 
repeated-measurement resting-state functional MRI 
to explore intersubject variability in connectivity. 
Individual differences in functional connectivity 
were heterogeneous across the cortex, with signifi- 
cantly higher variability in heteromodal association 
cortex and lower variability in unimodal cortices. In- 
tersubject variability in connectivity was significantly 
correlated with the degree of evolutionary cortical 
expansion, suggesting a potential evolutionary root 
of functional variability. The connectivity variability 
was also related to variability in sulcal depth but 
not cortical thickness, positively correlated with the 
degree of long-range connectivity but negatively 
correlated with local connectivity. A meta-analysis 
further revealed that regions predicting individual 
differences in cognitive domains are predominantly 
located in regions of high connectivity variability. 
Our findings have potential implications for under- 
standing brain evolution and development, guiding 
intervention, and interpreting statistical maps in 
neuroimaging. 


INTRODUCTION 


The human brain is characterized by striking interindividual vari- 
ability in neuroanatomy and function (Frost and Goebel, 2012; 
Rademacher et al., 2001; Sugiura et al., 2007; van Essen and 
Dierker, 2007) that is reflected in great individual differences in 
human cognition and behavior. Such variability is a joint output 
of genetic and environmental influences that may differentially 
impact on different brain systems (Glendenning and Masterton, 
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1998). For example, structural variability of association cortex 
is less influenced by genetic factors during development 
(Brun et al., 2009), allowing more variable impact of postnatal 
environmental factors that lead to the diversity of neural connec- 
tions beyond their genetic determination (Petanjek et al., 2011). 
A plethora of evidences suggest that neural systems subserv- 
ing higher-order association and integration processes are 
more variable than those implicated in unimodal processing. 
Language areas for example exhibit overproportionally high 
variability in cytoarchitectonically defined volume (Amunts 
et al., 1999), as well as in fMRI-derived localization (Frost and 
Goebel, 2012). At a macroscopic scale, structural variability in 
cortical folding is higher in association areas than in the motor 
cortex (Hill et al., 2010a). In addition, long association white 
matter fiber tracts are more variable than the optic radiation 
and the corticospinal tract (Burgel et al., 2006). In contrast to 
the large amount of work assessing structural variability across 
brain areas, individual variability in functional connectivity has 
not been systematically investigated and quantified. 

An individual brain might be best characterized by its connec- 
tome (Seung, 2012). One powerful technique for assessing 
connectivity utilizes fMRI data obtained under resting conditions, 
often referred to as intrinsic functional connectivity (Fox and 
Raichle, 2007). Individual differences in intrinsic functional 
connectivity can predict individual performance variability in 
several cognitive domains in the healthy (Andrews-Hanna 
et al., 2007; Seeley et al., 2007; van den Heuvel et al., 2009) 
and symptom severity in neuropsychiatric disorders (Fox and 
Greicius, 2010; Greicius, 2008). Quantifying the spatial distribu- 
tion of intersubject variability in connectivity could therefore 
provide new insights into the neural underpinnings of individual 
differences in human functions. This distribution could also 
have practical implications in guiding surgical mapping, inter- 
preting imaging results (if results are averaged across subjects, 
it is less likely to obtain a significant effect in highly variable 
regions) and understanding which areas are the most likely to 
relate to variability in behavior. 

In the present article, we collected intrinsic functional connec- 
tivity MRI data on 23 healthy subjects each scanned five times 
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Figure 1. Intersubject Variability in Resting-State Functional 
Connectivity Is Heterogeneous across the Human Cortex 
Intersubject variability was quantified at each surface vertex across 23 
subjects after correction for underlying intrasubject variability. Values below 
the global mean are shown in cool colors while values above the global mean 
are shown in warm colors. See also Figure S1. 


over 6 months. This unique data set allows us to assess the 
spatial distribution of intersubject variability while controlling 
for measurement instability based on intrasubject variance. 
This map of intersubject variability was then directly compared 
to maps of evolutionary cortical expansion, anatomical vari- 
ability, and long-range integration and regional segregation 
(Sepulcre et al., 2010). Finally we performed a meta-analysis to 
explore how functional connectivity variability may relate to 
previously observed individual differences in cognition and 
behavior. 


RESULTS 


Intersubject Connectivity Variability Ils Nonuniformly 
Distributed across Brain Networks 

Intersubject variability in intrinsic functional connectivity was 
quantified at each vertex of the brain surface after correction 
for nuisance variance (see Figures S1A and S1B, available on- 
line, and Experimental Procedures for the details). Intersubject 
variability demonstrated a nonuniform distribution across brain 
regions (Figure 1). Individual differences were largest in hetero- 
modal association cortex including the lateral prefrontal lobe 
and the temporal-parietal junction and minimal in unimodal 
sensory and motor cortices. Functional variability was also 
assessed within 7 specific brain networks (Yeo et al., 2011; 
Figure 2, top row). Intersubject variability within the boundary 
of each network was averaged and compared (Figure 2). We 
found that frontoparietal control and attentional networks 
demonstrated a high level of functional variability, whereas 
sensory-motor and visual systems were least variable. The 
default network demonstrated a moderate level of variability, 
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Figure 2. Functional Connectivity Variability Quantified across 
Cortical Networks 

The analysis was based on our prior parcellation of the cerebrum (Yeo et al., 
2011) into seven functional networks (top row), namely the frontoparietal 
control (FPN), ventral and dorsal attention (VATN, dATN), default (DN), limbic 
(LMB), sensory-motor (Mot), and visual (Vis) networks. Intersubject variability 
within the boundary of each network (black curves in the middle row) was 
averaged and plotted (bars in the bottom row). The dotted line indicates the 
global mean of intersubject variability in the entire cerebral cortex. See also 
Figure S2. 


which is lower than that of frontoparietal and attentional 
networks, but higher than the variability of sensorimotor and 
visual networks. 


Functional Connectivity Variability Ils Highly Correlated 
with Evolutionary Cortical Surface Expansion 

Functional connectivity variability was found to be highest in 
frontal, temporal, and parietal association cortex areas. These 
brain regions are phylogenetically late-developing regions 
(Kaas, 2006; Smaers et al., 2011) that are essential to complex 
and human specific cognitive functions like reasoning and 
language (Goldman-Rakic, 1988). As evolutionary history is 
usually represented by the phylogenetic tree, the fact that higher 
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Figure 3. Functional Connectivity Variability and Evolutionary 
Cortical Expansion Are Highly Correlated 

(A) The regional evolutionary cortical expansion between an adult macaque 
and the average human adult PALS-B12 atlas. Data were provided by van 
Essen and colleagues (van Essen and Dierker, 2007). On a whole-surface level, 
evolutionary expansion and functional variability (B) were significantly asso- 
ciated (r = 0.52, p < 0.0001). The correlation was shown in the scatter plot (C) 
where each 100" vertex is represented by a small circle. 


variability exists in the phylogenetically late regions may indicate 
an evolutionary root of the variability in functional connectivity. 
To test this hypothesis, we compared the functional variability 
map (Figures 1 and 3B) to a map of regional evolutionary cortical 
expansion between an adult macaque and the average human 
adult PALS-B12 atlas (Figure 3A) provided by David van Essen 
and colleagues (Hill et al., 2010b; van Essen and Dierker, 2007; 
http://sumsdb.wustl.edu/sums/directory.do?id=7601585). On a 
whole-surface level, evolutionary expansion and functional vari- 
ability were significantly correlated (r = 0.52, p < 0.0001; Fig- 
ure 3C), indicating that the extent of functional variability is 
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related to the evolutionary cortical expansion. To verify that 
this correlation is not influenced by the spatial dependence 
between neighboring vertices, we randomly sampled 7% of 
the vertices 1,000 times and computed the correlation coeffi- 
cient based on these subsets of vertices. These vertices were 
spatially independent as confirmed by Durbin-Watson test (See 
Experimental Procedures). All of the reported correlation coeffi- 
cients in our paper have been tested using this procedure and 
were not affected by spatial dependence between neighboring 
vertices. 


Functional Connectivity Variability ls Associated with 
Brain Folding Pattern but Not Cortical Thickness 

It has been well recognized that across individuals the cortical 
folding patterns are consistent in some regions but highly vari- 
able in some other regions (Hill et al., 2010a). Here, we investi- 
gated how the functional connectivity variability may relate to 
the known anatomical variability. Sulcal depth (see Experimental 
Procedures for definition and caveats) and cortical thickness 
were estimated for each subject using FreeSurfer (Figures 4A 
and 4B). To properly model the anatomical variability, we em- 
ployed the intraclass correlation (ICC; see Experimental Proce- 
dures) with the intrasubject variance sufficiently accounted for. 
Consistent with previous findings (Hill et al., 2010a), sulcal depth 
variability was most pronounced in lateral frontal and temporo- 
parietal regions but was low in the motor cortex. The default 
network showed moderate sulcal depth variability. In contrast, 
cortical thickness demonstrated a very distinct pattern with 
high variability in the motor area but low variability in the fronto- 
parietal network (See also Figure S3 for a quantification across 
seven functional networks). When quantified on the whole brain 
surface, sulcal depth variability showed a moderate but signifi- 
cant correlation with functional variability (r = 0.30, p < 0.0001), 
while cortical thickness variability was uncorrelated with func- 
tional variability (r = 0.05, p > 0.05). 


Functional Variability Is Positively Associated with the 
Degree of Long-Range Connectivity but Negatively 
Associated with Local Connectivity 

It has been suggested that developmental reorganization of 
functional connectivity is characterized by a shift of functional 
connectivity hubs from sensory-motor cortex toward default 
(Fransson et al., 2011) and frontoparietal network areas (Power 
et al., 2010). In adults, functional connectivity is known to form 
preferentially local connections within sensory and motor 
cortical regions, while hubs of distant connections are located 
in phylogenetically and ontogenetically later multimodal associ- 
ation cortices (Sepulcre et al., 2010). Here, we explored whether 
this special network organization of the human brain is related to 
functional variability. The degree of distant and local functional 
connectivity was quantified at each voxel in the brain volume 
according to Sepulcre et al. (2010). Distant connectivity was 
defined as the connection (r > 0.25) between two regions with 
a distance larger than 25 mm. Local connectivity was define as 
the connection (r > 0.25) within 12 mm. The percentage of distant 
connectivity (Figure 5A) demonstrated a moderate but significant 
correlation with the functional variability (r = 0.32, p < 0.0001) 
across the entire cerebral cortex (Figure 5B). Within the regions 
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dominated by local connectivity (blue regions in Figure 5A), func- 
tional variability showed a negative correlation (r = —0.33, p < 
0.0001) to the degree of local connectivity (Figure 5C). The rela- 
tion between functional variability and the degree of connectivity 
is exemplified in the default network. It has been reported that 
the default network is a hybrid hub of both local and long-range 
cortical-cortical interactions (Sepulcre et al., 2010). In our data, 
we have observed a moderate level of functional variability in 
the default network, consistent with the notion that functional 
variability is associated with the degree of long-range connec- 
tivity but negatively correlated with the degree of local 
connectivity. 





Regions Predicting Individual Differences in Cognitive 
and Behavioral Domains Are Predominantly Located in 
Regions of High Functional Connectivity Variability 
Intrinsic functional connectivity has been shown to reflect indi- 
vidual performance variability in several cognitive domains in 
healthy individuals (Seeley et al., 2007; van den Heuvel et al., 
2009). To determine if these regions previously shown to relate 
to individual differences in performance overlap with the 
currently identified regions of high intersubject variability, we 
performed a PubMed-based search of studies that reported 
associations between functional connectivity measures and 
individual differences in cognitive or behavioral domains 
including personality traits, memory performance, anxiety, risk 
seeking behavior, response inhibition, intelligence, and visual 
perception (for inclusion criteria, see Experimental Procedures; 
for a list of included studies, see Table S1). A total of 15 studies, 
comprising 573 subjects and 139 foci were retrieved. Quantifica- 
tion was performed on the brain surface (Figure 6) and the results 
revealed that about 73% percent of the clusters overlap with 
regions of high functional variability. Regions of high variability 
were defined as regions displaying variability above the global 
mean and covered about 51% of the cortical surface. 


Ruling Out Potential Confounding Factors 

To rule out the possibility that the observed functional connec- 
tivity variability was dominated by intersubject differences in 
head motion during the scan sessions, we calculated the mean 
relative displacement (Van Dijk et al., 2012) for each session of 
each subject. We chose a subset of ten subjects that displayed 
higher intra- than intersubject variance in head motion and quan- 
tified intersubject variability in functional connectivity using the 
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Figure 4. Relationship between Functional 
and Anatomical Variability 

Functional connectivity variability is significantly 
associated with the variability in sulcal depth (A) 
but not the variability in cortical thickness (B). 
Intersubject anatomical variability was calculated 
using intraclass correlation (ICC), with the intra- 
subject variance properly accounted for. Sulcal 
depth variability showed a significant correlation 
with functional variability (r = 0.30, p < 0.0001) 
while cortical thickness variability was uncorre- 
lated with functional variability (r = 0.05, p > 0.05). 
See also Figure S3. 


same procedure as in Figure 1. The functional variability map 
derived from this subset of subjects displayed the same charac- 
teristic topography as shown in Figure 1 (r = 0.77, p < 0.0001), 
suggesting that the functional variability obServed was not due 
to the intersubject variance in head motion. 

Higher degree of convolution in association cortex areas may 
also lead to lower fidelity of intersubject alignment in these 
regions (van Essen, 2005). To investigate this potential confound 
we regressed out sulcal depth variability, which comprises vari- 
ability due to alignment error, from the functional variability map. 
Figure S3 demonstrates that the overall pattern of functional 
connectivity variability remains stable after regression. Never- 
theless, this approach only partially accounts for alignment 
errors as it disregards cytoarchitectonic information of cortical 
areas, whose positions in relation to gyral and sulcal folds are 
themselves variable. We therefore further quantified functional 
connectivity variability in several histologically defined architec- 
tonic brain areas (Fischl et al., 2008) that are Known to show 
different susceptibility to misalignment. Previous studies have 
suggested that MT had a larger alignment error than BA 44/45 
(Yeo et al., 2010). We found that MT, although more prone to 
alignment variability, showed lower variability (0.60) in functional 
connectivity than BA 44/45 (0.64 and 0.65, respectively). This 
discrepancy may suggest that functional variability is influenced 
but not dominated by alignment variability. However, future 
investigation on architectonic variability across the brain will be 
useful to better address this potential confound. 

Two left-handed subjects had been included in our data set in 
order to roughly represent the handedness distribution in the 
healthy population. To investigate the potential impact of this 
handedness variability on our results, intersubject functional 
connectivity variability was recalculated after excluding the 
left-handed subjects (Figure S1C). The variability maps derived 
from these two data sets were highly correlated (r = 0.99), sug- 
gesting that the observed variability distribution was not domi- 
nated by the handedness variability in the data set. 


DISCUSSION 


Several findings in the current article add to our understanding of 
individual differences in functional connectivity. We demon- 
strated that functional connectivity variability has a specific 
topographic distribution with heteromodal association cortex 
being most variable and unimodal sensorimotor regions being 
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Figure 5. Functional Connectivity Variability ls Positively Associated with the Degree of Long-Range Functional Connectivity but Negatively 


Correlated with Local Connectivity 


Distant connectivity was defined as the connection (r > 0.25) between two regions with a distance larger than 25 mm. Local connectivity was defined as the 
connection (r > 0.25) within 12 mm. The percentage of distant connectivity was projected to the brain surface (A). Regions above the global mean are shown in 
yellow; regions below the global mean are shown in blue. Significant correlation (r = 0.32, p < 0.0001) was found between the functional connectivity variability and 
the percentage of distant connectivity across the entire cerebral cortex (B). In the regions dominated by local connectivity, functional connectivity variability was 
negatively correlated (r = —0.33, p < 0.0001) with the degree of local connectivity (C). 


least variable. This functional connectivity variability is related to 
evolutionary cortical expansion and variability in cortical folding 
pattern but not cortical thickness. Further analyses revealed 
that functional connectivity variability is associated with network 
properties of functional integration and segregation. Finally, we 
demonstrated that our map of functional connectivity variability 
overlaps well with prior reports linking individual differences in 
functional connectivity to behavioral performance. 


Potential Causes of Strong Variability in Association 
Cortex 

Functional variability in human cerebral cortex is likely to be the 
result of evolution that has shaped a unique distribution of 
susceptibility to genetic and environmental influences. Associa- 
tion cortex areas, where functional architecture appeared to be 
most variable, are phylogenetically late-developing regions that 
underwent a disproportionate enlargement during human evolu- 
tion (Kaas, 2006; Smaers et al., 2011; van Essen and Dierker, 
2007). Evolutionary trajectories can be partially retraced in indi- 
vidual development (Clancy et al., 2000), where association 
cortex exhibits the most protracted course of white (Yakovlev 
and Lecours, 1967) and gray (Gogtay et al., 2004; Shaw et al., 
2008) matter maturation and most pronounced postnatal cortical 
expansion (Hill et al., 2010b). This prolonged maturation course 
of association cortex implicates a prolonged exposure to vari- 
able extrinsic experience during a time of high neuroplasticity 
(Petanjek et al., 2011). In addition, structural variability of late 
maturing association cortex is probably less genetically influ- 
enced during development (Brun et al., 2009), again enabling 
more variable impact of postnatal environmental factors that 
lead to the diversity of neural connections beyond their genetic 
determination (Petanjek et al., 2011; but also see Chen et al., 
2011; Rimol et al., 2010; Thompson et al., 2001 for different 
explanations). Besides this prior evidence on heritability of 
anatomical properties, functional connectivity, e.g., of the 
default network, is known to be influenced by genetic factors, 
which cannot necessarily be attributed to anatomical variability 
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(Glahn et al., 2010). Nevertheless, the spatial distribution of the 
heritability of functional connection strength across the entire 
brain is yet to be unveiled. Finally, the dynamics of synaptic over- 
production in early childhood and consecutive synaptic pruning 
may contribute to a similar functional hierarchy, where synaptic 
overproduction is highest in the prefrontal cortex and lowest in 
primary sensory regions (Elston et al., 2009; Jacobs et al., 
1997). High synaptic overproduction may provide more freedom 
for selective stabilization to operate on during development. 
Taken together, a protracted maturation, weaker genetic influ- 
ence on structure and more synaptic over-production may jointly 
contribute to the high functional variability of multimodal associ- 
ation cortices as reported in this study. 


Functional Variability Ils Related to Evolutionary Cortical 
Expansion and Cortical Folding 

Functional variability is correlated with variability of sulcal depth, 
a proxy of cortical folding. From an evolutionary perspective, the 
degree of gyrification is highest in phylogenetically young asso- 
ciation cortex and lowest in phylogenetically older occipital and 
motor cortex (Zilles et al., 1997), resulting in highest sulcal depth 
variability and positional variability in association cortex (Hill 
et al., 2010a). It is striking how well the regional evolutionary 
cortical expansion and sulcal depth variability maps match the 
distribution of variability in functional connectivity as revealed 
in this study. In contrast, no significant correlation was found 
between cortical thickness variability and functional variability. 
This finding is consistent with the fact that brain evolution has 
been characterized by huge surface expansion (e.g., ten-fold 
between macaque and human [Preuss, 1995]) without a signifi- 
cant increase in cortical thickness (Rakic, 1995). 


Functional Variability Is Related to the Need for 
Long-Range Information Exchange 

The human brain possesses a complex architecture with some 
areas highly specialized for local, modular processing and 
certain areas connecting and integrating these otherwise 
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Figure 6. Loci that Predict Individual Differences in Behavioral and 
Cognitive Domains Are Predominantly Located in Cortical Areas of 
High Functional Connectivity Variability 

Loci were derived from a meta-analysis that included 15 studies that found 
associations between functional connectivity and individual differences in 
cognitive and behavioral domains. Loci were merged in the volume. Frequency 
of contributing foci was estimated for each voxel. Results were smoothed, 
normalized, and projected to the surface. Quantification revealed that about 
73% of the clusters associated with individual differences are located in 
cortical regions that display high functional connectivity variability (above the 
global mean, 51% of the cortical surface in total). For an overview of included 
studies, see Table $1. 


segregated brain regions or systems (Buckner et al., 2009; 
Sepulcre et al., 2010). Such arrangement may maintain high 
information processing efficiency given that the human brain 
has tripled the size over past several million years. The ratio of 
local to distributed areal projections is suggested to be critical 
to the evolution of higher-order cognitive functions including 
language, reasoning, and foresight (Kaas, 2005; Semendeferi 
et al., 2001). We have previously reported that regions within or 
near primary sensory and motor areas display high local connec- 
tivity consistent with a modular organization. In contrast, distant 
connectivity is prominent across association areas in parietal, 
lateral temporal, and frontal cortices as well as paralimbic cortex 
including posterior cingulate (Sepulcre et al., 2010). Here, we 
extend these insights by showing that functional variability is 
strongly correlated with the degree of distant connectivity but 
negatively correlated with the degree of local connectivity. A 
potential inference from this observation is that functional vari- 
ability may not become prominent until distant connectivity 
emerges, i.e., species with smaller brains dominated by local, 
modular processing may have limited functional variability, 
hence the more uniform and predictable behavior. 

A particularly intriguing observation is that the default network, 
which represents a hybrid hub for local and distant connections 
(Sepulcre et al., 2010), exhibited moderate variability in functional 
connectivity across subjects. Given that the DMN has been 





described as a network that subserves mostly internal thought 
processes and human specific functions such as autobiograph- 
ical memory retrieval (Svoboda et al., 2006), imagining the future 
(Schacter et al., 2007), and mind wandering (Mason et al., 2007) 
one would predict the DMN to be among the most variable 
networks, both across and within subjects. Neither of those two 
predictions was proven correct in this study. The DMN showed 
low intrasubject variability (Figure S1B) and intermediate inter- 
subject variability (Figures 2 and S2). Taking into account that 
the DMN is present in rodents (Lu et al., 2012) and anaesthetized 
monkeys (Vincent et al., 2007) it seems plausible that the DMN 
subserves both phylogenetically older, putatively less complex 
functions and human specific higher order cognitive functions. 
This could be reflected in the intermediate variability of the DMN 
where information processing that involves modular computa- 
tion could be consistent across subjects, whereas processing 
which associates distributed information from key limbic, parietal, 
and prefrontal regions exhibits strong individual differences. 


Clinical Relevance of Individual Differences 

While functional variability in association cortex has important 
implications for the evolution of higher-order cognitive abilities, 
it might also relate to an increased susceptibility to the formation 
of abnormal circuitry as manifested in neuropsychiatric disor- 
ders. Here, we demonstrate that individual differences in mental 
domains such as personality traits can be linked to brain regions 
of high functional variability. A caveat is that based on the avail- 
able literature the majority of included studies investigated indi- 
vidual differences in higher cognitive functions, which might 
constitute a publication bias favoring higher order association 
cortex areas in displaying individual differences. 

Functional development of the human brain is characterized 
by a general trend toward increases in connectivity across 
widely distributed regions, conceptualized as the development 
of a ‘local to distributed’ organization (Fair et al., 2009). Studies 
have suggested that abnormal development leading to variable 
disconnection of focal brain regions, especially regions that are 
functionally integrated network hubs, might be present in many 
neuropsychiatric disorders (Zhang and Raichle, 2010). In this 
context, it is noteworthy that many neuropsychiatric diseases 
such as anxiety disorders, bipolar disorder, depression, eating 
disorder, psychosis (including schizophrenia), and substance 
abuse most commonly emerge during adolescence (Kessler 
et al., 2007), a period critical for the establishment of long-range 
connection hubs that signify functional variability. Brain circuits 
susceptible to neuropsychiatric diseases may therefore be iden- 
tified based on the abnormal range of connectivity variability in 
patients. Knowing cortical areas of highest individual variability 
may furthermore help guide investigations into individual differ- 
ences in disease susceptibility. 

As clinical practice moves ever closer to the goal of individual- 
ized therapy, Knowing the distribution of individual differences in 
brain connectivity is likely to be important. For example, func- 
tional and connectivity data are playing an increasing role in 
guiding operative approaches (Liu et al., 2009). Knowing that 
a surgical resection is near an area of high intersubject variability 
may inform acquisition of preoperative imaging. Similarly, there 
is increasing evidence that therapeutic brain stimulation might 
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be guided by differences in connectivity (Fox et al., 2012). 
Knowing whether a target of brain stimulation is in an area of 
high or low individual variability will be important for determining 
whether one can target based on group averages or if one should 
obtain information on a patient’s specific connectivity pattern. 


Relevance to Interpretation of Group Maps in 
Neuroimaging 

Regardless of whether one is studying functional connectivity 
or task-based activations, neuroimaging results are generally 
presented as a statistical map computed across a group of 
subjects. The creation of these statistical maps necessarily 
incorporates the variance across the group. As such, the map 
of individual differences presented here is highly relevant for 
interpretation of these statistical images. Specifically, one is 
more likely to get a significant result in areas of low individual 
variability such as primary sensory or motor cortex and less likely 
to get a significant result in areas of high individual variability. 
Therefore, the risk of false-positives and false-negatives in neu- 
roimaging is likely non-uniformly distributed across the human 
cortex. Variance maps from an independent data set such as 
the one presented here might eventually be used to formally 
correct for this heterogeneity in creating statistical images. 


EXPERIMENTAL PROCEDURES 


Participants and Data Collection 

Twenty-five healthy subjects (age 51.8 + 6.99, 9 female) were recruited for 
a longitudinal feMRI study. The data was collected as a control sample of 
a longitudinal stroke study. Therefore the age range is slightly higher than 
what would be expected for a study of healthy adult subjects. The data 
set also included two left-handed subjects, roughly representing the handed- 
ness distribution in the healthy population (Connolly and Bishop, 1992). 
Participants were screened to exclude individuals with a history of neurologic 
or psychiatric conditions as well as those using psychoactive medications. 
Participants provided written informed consent in accordance with guidelines 
set by institutional review boards of Xuanwu Hospital. Each subject under- 
went five scanning sessions within 6 months (7, 14, 30, 90, and 180 days 
from the enrollment). All participants performed two or three rest runs per 
session (6 m 12 s per run) to estimate intrinsic functional connectivity. After 
quality control, 23 subjects who had at least two good runs (tSNR > 100) in 
each session were included in this study (mean = 2.02 runs). All data were 
acquired on a 3 Tesla TimTrio system (Siemens) using the 12-channel 
phased-array coil supplied by the vendor. Functional data were obtained 
using a gradient echo-planar pulse sequence (TR, 3,000 ms; TE, 30 ms; flip 
angle, 90°; 3 mm isotropic voxels, transverse orientation, 47 slices fully 
covering cerebral cortex and cerebellum). Subjects were instructed to stay 
awake and keep their eyes open; no other task instruction was provided. 
Structural images were acquired using a sagittal MP-RAGE three-dimensional 
T1-weighted sequence (TR, 1600 ms; TE, 2.15 ms; flip angle, 9°; 1.0 mm 
isotropic voxels; FOV, 256 x 256). 


Data Preprocessing 

Resting-state fMRI data were processed using previously described proce- 
dures (Van Dijk et al., 2010; Yeo et al., 2011). Structural data was processed 
using the FreeSurfer version 4.5.0 software package (http://surfer.nmr.mgh. 
harvard.edu). Surface mesh representations of the cortex from each individual 
subject’s structural images were reconstructed and registered to a common 
spherical coordinate system (Fischl et al., 1999). The structural and functional 
images were aligned using boundary-based registration (Greve and Fisch, 
2009). The resting-state BOLD fMRI data were then aligned to the common 
spherical coordinate system via sampling from the middle of the cortical ribbon 
in a single interpolation step. See Yeo et al. (2011) for details. 
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In this study, a symmetric surface template of the cerebral cortex (unpub- 
lished) was constructed using FreeSurfer. fMRI data of each individual were 
then registered to this template. The data were resampled on this template 
with a mesh of 1,284 vertices. For each vertex in this mesh, the nearest vertex 
in the higher resolution template was extracted and if multiple nearest vertex 
existed, the values on these vertices were averaged. We have used this lower 
resolution template to achieve computational efficiency but this re-sampling 
procedure may introduce noise. This was mitigated by a smoothing prepro- 
cessing step that we have taken. 


Estimating Intersubject Functional Variability 
Functional correlation maps were computed by taking each of the 1,284 
vertices as the seed, resulting in 1,284 maps for each subject and session. 
The correlation map based on each seed vertex can be denoted as F;(s,t), 
where /=1,2,...1284, and F; is a1 X 1284 vector, s indicates the subject, ft 
indicates the session. 

For a given seed vertex /, the similarity between the 23 maps derived from 
23 subjects was quantified by averaging the correlation values between any 
two maps: 


Ri(t) =E|corr(Fi (Sp, t), Fi(Sq,t))], where p,q =1,2, ...23;p#q. 


The intrasubject variance was estimated using the 5 maps derived from 5 
scanning sessions of each subject: 


N,(s) =1 — E[corr(F;(s, tm), Fi(s,t,))], where m,n=1,2,...5;m#n. 


The intrasubject variance was then averaged across 23 subjects and 
assigned to the seed vertex / (see Figure S2 top row): 


N; = E(N; (s)| F 


Note that the intrasubject variance consists of the variance caused by 
technical noise, which may be reflected by the tSNR of the BOLD signal 
(Figure S1B, middle row), as well as the biological variance related to the brain 
state change within subjects (Figure S1B, bottom row). 

To estimate intersubject variability, the similarity map A;(t) was first inverted 
(by subtraction from 1; see Figure S1A) and then the intrasubject variance was 
regressed out using ordinary least-squares regression (i.e., a general linear 
model, GLM). The residual map was taken as the estimate of functional 
variability, 


Vi(t) =[1 — Ri(t)] — BN; —c, 


where 6 and c are parameters determined via ordinary least-squares. Vari- 
ability maps derived from each session t are averaged and shown in Figure 1. 


Parcellation and Seed-Based Network Analysis 
To quantify variability in specific functional networks, we used the functional 
atlas derived from a clustering approach (http://surfer.nmr.mgh.harvard.edu/ 
fswiki/CorticalParcellation_Yeo2011; Yeo et al., 2011). The boundaries of 
seven networks were projected to the symmetric surface template. Intersub- 
ject variability values were then averaged within each network (Figure 2A). 
For the ROIl-based analysis described in Figure S2, we used a group of 
regions of interest (ROI) associated with different brain functions (Van Dijk 
et al., 2012), including the ATN (paraCG, FEF, MFG, Insula, MTplus, TPu, 
sPL), the FPN (aPFC, ApCC, dIPFC, SFG, IPL), the DN (PCC, aMPFC, dMPFC, 
LPC, SFG, LTC), the motor cortex (hand, foot, tongue region), the visual cortex 
(medial and lateral V1), and the auditory cortex (STG). Seeds were created by 
projecting the center of each volume ROI! (MNI152 volumetric space) to the 
FreeSurfer spherical surface model and constructing a circle (radius = 8 mm, 
defined as the arc length on the sphere) around each projected peak vertex 
on the sphere. The coefficient of variance of the correlation strength between 
a given pair of seeds was computed as the standard deviation divided by the 
mean across 23 subjects. To account for the measurement instability, the raw 
intersubject variance was normalized by the mean intrasubject variance for 
each seed pair. The coefficient of variance was then averaged across all 
seed pairs of one network. The surface-based ROIs may correspond to 
different sizes of brain volume but this source of variability is not significantly 
affecting the result. The ROl-based analyses described above were repeated 
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using standard volumetric spherical seeds in the volumetric space, the re- 
ported ranking of variability among functional networks remained unchanged. 


Relation to Evolutionary Cortical Expansion 

The map of regional evolutionary cortical expansion between an adult 
macaque and the average human adult PALS-B12 atlas was published previ- 
ously (Hill et al., 2010b; van Essen and Dierker, 2007) and made publicly avail- 
able. The right hemisphere evolutionary expansion map and the functional 
variability map were projected to the Conte69 164k_fs_LR mesh (van Essen 
et al., 2012) (http://sumsdb.wustl.edu/sums/directory.do?id=8291494&dir_ 
name=CONTE69). The data were extracted using the Caret Surface Statistics 
Toolbox (Diedrichsen, 2005) for the correlation analysis. The absolute expan- 
sion ratio was normalized by taking the logarithm and subtracted with a 
constant. 


Relation to Anatomical Variability 

Sulcal depth and cortical thickness measurements were calculated using 
FreeSurfer (Fischl et al., 1999). The sulcal depth estimated by FreeSurfer is 
not a direct measure of distance to the outer cortical margin, but the integrated 
dot product of the movement vector with the surface normal during inflation. It 
highlights large-scale geometry as deep regions consistently move outward 
and get a positive value while superficial regions move inward and get a nega- 
tive value. Intersubject variability in sulcal depth and cortical thickness was 
estimated vertex-wise using intraclass correlation (Shrout and Fleiss, 1979) 
with the intrasubject variance accounted for. The Pearson’s correlation 
coefficient was calculated between functional variability and anatomical vari- 
ability across the whole brain. To demonstrate the topological impact of 
anatomical variability on functional variability, a GLM approach was applied 
to regress out sulcal depth and cortical thickness variability from the functional 
variability map. 


Testing the Potential Impact of Spatial Dependence on Correlation 
Analyses 

To test the potential impact of spatial dependence between neighboring 
vertices on correlation analysis, we performed a repeated (n = 1,000) random 
sampling of 7% of the vertices and computed the correlation coefficient on the 
subsets of the vertices. For each subset, the Durbin-Watson test was per- 
formed to estimate the spatial dependence (DW > 2). Correlation coefficients 
were averaged across the 1,000 iterations. 


Meta-analysis of Individual Differences Predicted by Functional 
Connectivity 

We performed a voxel-wise frequency-based meta-analysis. A PubMed 
search was conducted using three sets of search terms: (1) search: individual 
differences, intrinsic connectivity; (2) search: individual differences, resting- 
state fMRI; (8) search: individual differences, connectivity, MRI. After 
accounting for redundancies, this resulted in 182 studies to be reviewed. 
The following inclusion criteria were applied: healthy, adult human subjects, 
original research, fMRI study, reported cerebral/cortical coordinates in stan- 
dardized stereotaxic space (Talairach or Montreal Neurological Institute 
[MNI] template) and association between an individual cognitive/behavioral/ 
psychological trait and a functional connectivity measure. Fifteen studies 
met inclusion criteria. The meta-analysis was conducted in MNI space. For 
studies that reported coordinates in Talairach space (Talairach and Tournoux, 
1988), used SPM or FSL, and did not specify the use of Lancaster transforma- 
tion (Laird et al., 2010), conversion to MNI coordinates was performed using 
the reversed Brett transformation (Brett et al., 2001). For studies that reported 
coordinates in Talairach space and used neither SPM nor FSL, conversion to 
MNI coordinates was performed using the (reversed) Lancaster transformation 
(Laird et al., 2010). Three millimeter soheres around each focus were merged 
in MNI 152 volumetric space. For each voxel, the number of contributing foci 
was calculated. The resulting volume map was spatially smoothed (FWHM 
12 mm), normalized (z score) and projected to the surface. 


Visualization 
Since the main analyses were performed in FreeSurfer symmetric surface 
space, the final results of both hemispheres were projected only to the left 





hemisphere of the inflated PALS cortical surface using CARET (van Essen, 
2005) for the purpose of visualization. The right hemisphere results shown in 
the figures were mirrored from the results rendered on left CARET surface. 
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